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Abstract. We consider several aspects of non-periodic Ising models in one 
and two dimensions. Here we are not interested in random systems, but 
rather in models with intrinsic long-range aperiodic order. The most promi- 
nent examples in one dimension are sequences generated by substitution 
rules on a finite alphabet. The classical one-dimensional Ising chain as well 
as the Ising quantum chain with coupling constants modulated according 
to substitution sequences are considered. Both the distribution of Lee-Yang 
zeros on the unit circle in the classical case and that of fermion frequen- 
cies in the quantum model show characteristic gap structures which follow 
the gap labeling theorem of Bellissard. We also investigate the zero-field 
Ising model on two-dimensional aperiodic graphs, which are constructed 
from rectangular grids in the same spirit as the so-called Labyrinth. Here, 
duality arguments and exact solutions in the sense of commuting transfer 
matrices are used to gain information about the critical behaviour. 



1. Introduction 

The study of lattice spin models has long been a major subject of statistical 
physics. These systems, of which the Ising model is among the simplest, are 
interesting both from a mathematical point of view and for applications in 
physics. In particular, the universality hypothesis states that the behaviour 
of such systems close to critical points (provided the interactions are of 
finite range) depends only on very general properties as, for instance, the 
dimension and symmetries of the model. This means that real systems in 
the vicinity of second-order phase transitions (for example magnets at the 
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Curie temperature) may well be described by very simple toy models which 
only capture the essential ingredients. 

Unfortunately, even such a simple model as the Ising model has not 
been solved in three dimensions, which clearly is the physically most in- 
teresting case. However, ID and 2D models not only initiated important 
new developments in mathematics and theoretical physics, but also found 
various physical applications (Stinchcombe, 1973; Einstein, 1987). 

While there is a wealth of detail and rigorous treatment of periodic 
and related systems (Ruelle, 1969; Ellis, 1985; Simon, 1993), the knowledge 
about random systems is still rather vague and that about aperiodic sys- 
tems with long-range orientational order hardly existing. One exception is 
the article by Geerse and Hof (1991), where an approach to equilibrium 
states for lattice gas models on self-similar tilings is made, but, apart from 
that, few exact results about the periodic case have found an aperiodic 
counterpart. 

In this article, we present results on one- and two-dimensional aperiodic 
Ising models, where the aperiodicity is implemented by means of sequences 
based on substitution rules. We should emphasize that we do not aim at 
giving final answers to the problems connected with aperiodicity in lattice 
models, but rather collect useful insight for further, and more rigorous, 
investigations. 

The article is organized as follows. After brief remarks about the his- 
tory of the Ising model. Section 2 gives a short review of the properties 
of (periodic) Ising models in one and two dimensions. At the same time, 
this part introduces our notation and some basic concepts. In Section 3 we 
introduce the classical ID Ising model (in constant magnetic field) with an 
aperiodic sequence of coupling constants. We investigate its Lee-Yang or 
magnetic field zeros which (in proper variable) lie on the unit circle but 
show some unexpected behaviour: generically they only fill a Cantor subset 
of it (a first step of a proof is given in the Appendix). 

Next, Section 4 deals with its quantum mechanical counterpart, the 
Ising quantum chain with transversal field. It corresponds to an anisotropic 
limit of the the classical 2D Ising model without magnetic field (Fradkin 
and Susskind, 1978; Kogut, 1979). Starting from the Harris (1974) criterion 
for relevant disorder and Luck's (1993) adaption to quantum chains, we 
discuss a selection of typical examples with respect to critical behaviour 
and conformal invariance. It turns out that the latter is lost if fluctuations 
in the couplings are relevant. 

In Section 5 we move on to classical 2D Ising models (in zero external 
field) on non-periodic graphs. By means of Baxter's (1978) concept of "Z- 
invariance" we derive an exact solution for a subclass of models with a phase 
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transition of the Onsager universality class. In contrast to the situation in 
Section 4, this result is robust to a larger class of fluctuations. 

Though this article is not a review, we have spent some effort to include 
a rather complete list of articles (not all of which are cited in the text) in 
the Bibliography at the end of this paper. 

2. The Ising Model 

The Ising model is, without doubt, the most famous and best understood 
among the models of statistical mechanics. It was proposed by Lenz (1920) 
as a toy model of ferromagnetism. Five years later, his student Ising (who re- 
cently celebrated his 95th birthday) published the solution of the ID model 
(Ising, 1925), which showed no phase transition at finite (i.e., non-zero) 
temperature, and Ising supposed the same to be true in any dimension. 

This conclusion, however, proved wrong. In fact, Peierls (1936) showed 
that in two and more dimensions the Ising model undergoes a phase tran- 
sition from a disordered high-temperature to an ordered low-temperature 
phase. By duality arguments, Kramers and Wannier (1941) located the 
critical point of the 2D zero-field square lattice Ising model, assuming that 
the phase transition is unique. Shortly after, Onsager (1944) published his 
famous solution for the partition function of this model. The spontaneous 
magnetization in the ordered phase and its power law behaviour with the 
critical exponent (3 = 1/8 was announced by Onsager and calculated soon 
after (Yang, 1952; Montroll et al., 1963). Correlation functions have also 
been studied, see for instance (McCoy and Wu, 1973). For further details 
on the history and an extensive bibliography we refer to (Brush, 1967). 

Today, inspired by Onsager's work (Baxter, 1995), a plethora of so- 
called "solvable models" in two dimensions is known (Jimbo, 1989). These 
are constructed as solutions of the Yang-Baxter equation (YBE) which ba- 
sically is Onsager's star-triangle relation (see (Au-Yang and Perk, 1989) 
for its history). The YBE implies that row transfer matrices (Baxter, 1982) 
form a one-parameter family of commuting matrices. In a sense, the corre- 
sponding models can be regarded as generalizations of the Ising model, but 
for most of them solutions are only known (under a list of assumptions) in 
the thermodynamic limit whereas Onsager was also able to compute the 
complete finite-size spectrum of the Ising model. 

To the present day, the 2D Ising model in a magnetic field has not been 
solved. However, some properties of this system are known from conformal 
field theory (Zamolodchikov, 1989). Moreover, a solvable model has recently 
been found (Warnaar et al., 1992) that (in the sense of universality) de- 
scribes the same physics as the Ising model at the critical temperature in 
the presence of a symmetry-breaking field. 
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By introducing disorder (for instance by considering Ising models with 
random coupling constants) into a model one might change its physical 
properties considerably, see (Koukiou, 1993; McCoy, 1995) and references 
therein. There are general qualitative arguments about the effects of disor- 
der. Clearly, one would expect that critical singularities, if affected at all, 
can only be "softened" by randomness because these are co-operative phe- 
nomena. This is precisely what was observed by McCoy and Wu (1968): a 
certain kind of layered disorder destroys the logarithmic singularity of the 
specific heat at the phase transition of the 2D Ising model. A well-known 
criterion for the relevance of disorder to critical properties is due to Harris 
(1974) which depends on the sign of the specific heat exponent a in the un- 
perturbed model: critical singularities are smoothed out for a > 0, whereas 
disorder is irrelevant for a < 0. However, since a logarithmic divergence 
of the specific heat means a = 0, this criterion does not give a definite 
answer for the 2D Ising model. Numerous investigations of random bond 
Ising models using analytical and numerical techniques have not come up 
with a definite answer either, but in many cases no differences between the 
random and periodic case could be detected. 

In what follows, we are interested in a kind of aperiodicity which is dif- 
ferent from randomness, and in fact constitutes a long-range order in the 
system (wherefore the term "disorder" is not appropriate). Typical exam- 
ples are Ising models with coupling constants chosen according to (com- 
pletely deterministic) substitution sequences, or Ising models on quasiperi- 
odic graphs such as the Penrose tiling. In particular for the latter there is 
an extensive literature of mainly numerical results, mostly based on Monte 
Carlo simulations. To make some progress with open questions, we focus on 
exact results. We now set the scene by briefly reviewing the periodic case. 



2.1. ID ISING CHAIN 

Consider a linear array of N "Ising spin" variables aj, enumerated by j = 
1, . . . , iV (see Figure 1), which take values G {1, -1} (or "-|-" and "-", for 

short). To any configuration a = {uj | j = 1, • • • , G {1, —1}^ a quantity 
[energy) is assigned through 

N N 

^w = -E'^.^.^.+i-E^.^. (1) 
j=i j=i 

where it is convenient to use periodic boundary conditions, i.e., (Tjy_|_^ = 
(7^. For the moment, we keep our notation general and allow the coupling 
constants Jj and the local values of the magnetic field Hj to depend on 
the position in the chain. The model is called ferromagnetic if the energy 
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12 3 N-1 N 

Figure 1. One-dimensional Ising chain 

weight favours the variables on neighbouring sites to coincide, i.e., if the 
coupling constants Jj are positive. Although many results remain valid 
for the antiferromagnet as well, we will usually assume that all couplings 
are ferromagnetic. In particular, we do not discuss Ising models with a 
mixture of both ferro- and antiferromagnetic couplings which already in 
one dimension show complicated frustration effects, see e.g. (Luck, 1987). 
The statistical or Boltzmann weight of a configuration a is 

N 

W{a) = n W^{a^, cT^.+i) = exp(-/3^(a)) (2) 
j=i 

which is a product of local Boltzmann weights 

(d, a') = exp (^K^aa' + ^ {h^a + h^+i(j')^ (3) 

associated to the bond between sites j and j + 1. Here, we introduced 
dimensionless quantities Kj = (3Jj and hj = (3Hj, where (3 = l/k^T with 
temperature T and Boltzmann's constant kg. 

The (canonical) partition function is defined by 

Z^ = J2^xp{-/3E{a)) (4) 

where the sum is over all 2^ configurations a. Then W{a)/Zj^ defines 
a probability measure (a finite volume Gibbs state) and expectation val- 
ues {correlation functions) are computed with respect to this measure (for 
example the local magnetization {a-) = Y^^a-W{a)/Zj^ at site j). Intro- 
ducing local transfer matrices Tj defined by 

[w]{-,+) W]{-,-) I ^ ( e-^'^'j-(''j-''j + i)/2 ^Kj-{hj+hj+i)/2 j 

allows us to rewrite the sum over configurations in the partition function 
in terms of matrix multiplications 

N 

% = E n ^.(^.' = ^'(T,T, ■■■T^)= tr(r) (5) 

<T j = l 
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where the trace is due to periodic boundary conditions. The matrix T is 
frequently referred to as monodromy matrix in the literature. The (specific) 
free energy is given by 

/iV = -^log(^iv) (6) 

from which other thermodynamic quantities result as suitable derivatives. 

One is particularly interested in the thermodynamic limit in which the 
system size becomes infinite, and especially in points where the resulting 
free energy develops non-analytic behaviour as a function of temperature 
which can be regarded as a definition of a phase transition (Griffiths, 1972; 
Thompson, 1979). Clearly, if the local weights are analytic functions of 
temperature (as is usually the case), phase transitions can only appear 
in the thermodynamic limit. Of particular interest are so-called critical 
points (or second-order phase transitions) which are those phase transition 
points where the first derivative of the free energy is still continuous, while 
the second is not. These are characterized by correlation functions which 
show algebraic (with powers that are known as critical exponents or crit- 
ical indices) rather than the usual exponential decay (with characteristic 
distances called correlation lengths). In other words, correlation lengths di- 
verge if one approaches a critical point, which means that fluctuations (not 
just microscopically, but at all length scales) become important and crit- 
ical systems can effectively be described by continuum field theories, see 
Section 2.4. 

In general, it is impossible to calculate the partition function (5) explic- 
itly, even for this simple ID model. However, if all the transfer matrices 
Tj commute with each other, we can diagonalize them simultaneously. For 
example, for the translationally invariant case Jj = J and Hj = H, one 
obtains = T for all j (hence T = T^) and 

/. = -5^^1og(MT-)) = -5^^1og(Af + Ar) (7) 

where A^ > Aj denote the two eigenvalues of the (positive) symmetric 
matrix T . The free energy / = limjv->oo /jv the thermodynamic limit is 

-(3f = log(Ai) = log (^e^^ cosh(/i) + ^e2^^' sinh2(/i) + ^ (8) 

and hence is an analytic function of the magnetic field H and the tem- 
perature T for all real H and positive T. The point H = T = can be 
interpreted as sort of a critical point (in the sense that the correlation 
length becomes infinite), see e.g. (Baxter, 1982, pp. 32-38). The situation 
is generic for ID models with short-ranged interactions which cannot have 
phase transitions at non-zero temperature (van Hove, 1950). 



APERIODIC ISING MODELS 



7 




Figure 2. Magnetic field zeros of an alternating Ising chain and their integrated density 

Before we move on to 2D systems, let us have a short look at the parti- 
tion function Zjy (5) for an Ising chain with uniform magnetic field Hj = H 
and its zeros in the variable uj = exp{2(3H) {Lee-Yang zeros). According to 
the Lee-Yang theorem (Lee and Yang, 1952; Ruelle, 1969), for ferromag- 
netic coupling constants the zeros lie on the unit circle. For the simply 
periodic Ising chain {Jj = J), the zeros in the thermodynamic limit fill a 
connected part of the unit circle densely (though not uniformly) leaving a 
gap around the intersection with the positive real axis which only closes 
for r — 7- when the distribution of zeros becomes uniform. More generally, 
one obtains p connected parts of zeros on the unit circle for periodic chains 
with unit cell of length p. 

In Figure 2, we show the distribution of the magnetic field zeros on 
the unit circle and their integrated density (plotted against the angle in 
units of 27r) for an alternating Ising chain {p = 2) with coupling constants 
J2j-i = Ja ^^'^ '^2j — '^b^ with numerical values exp{2(3J^) = 3/2 and 
exp(2/3J^) = 3 (i.e., Jb/Ja ~ 2.71). The derivation is given in the Appendix. 

2.2. SQUARE LATTICE ISING MODEL 

Let us have a short look at the Ising model on the square lattice. In Figure 3, 
the Ising variables a^^, G {1,-1} are represented by black dots on the 
vertices of the square lattice, neighbouring pairs of spins interacting with 
coupling constants K/(3 along horizontal and L/(3 along vertical bonds, 
respectively. The energy for the Ising model in a magnetic field H = h/(3 
on a lattice of size NxM reads 

N M 

(3E{a) = -J2J2 {l^('j,k(^j+i,k + + (9) 

j=i k=i 
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Figure 3. The Ising model on the square lattice and its interpretation as an IRF model 



where periodic boundary conditions in both horizontal and vertical direc- 
tions are assumed (thus we actually work on a torus). 

For the solvable zero-field case, there is a variety of approaches to cal- 
culate the free energy, see (Baxter, 1982) and references therein. Let us 
add how the model fits into the zoo of solvable lattice models in the sense 
of commuting transfer matrices. Figure 3 also shows another square lat- 
tice of dashed lines which contains additional vertices represented by open 
circles. Using this lattice, the Ising model can be reformulated as an IRF 
{interaction-round- a- face) model, where statistical weights are associated 
to elementary plaquettes (faces) rather than to bonds of the lattice. To 
do this, one defines a three-state model (with local states {1,0, —1}, say) 
on the dashed lattice subject to the requirements that states on adjacent 
vertices of the lattice have to differ by one. In particular, the states 1 and 
— 1 may not be neighbours on adjacent lattice sites (or, equivalently, the 
statistical weights of such configurations are zero). Then any "allowed" con- 
figuration has one sublattice filled with the state 0, while only 1 and — 1 
occur on the other sublattice. One can now find face weights for this model 
such that its partition function coincides with that of the Ising model, see 
e.g. (Wadati et al., 1989). For the zero-field case H = 0, these face weights 
can be parametrized by elliptic functions of a so-called spectral parameter 
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such that they satisfy the YBE, which implies commutativity of transfer 
matrices for different spectral parameters. So, eigenvalues depend on the 
parameter but eigenvectors do not. 

The physical quantity (often called order parameter) that distinguishes 
between the ordered and disordered phases of the ferromagnetic Ising model 
is the magnetization 

. N M 

m(H,T)= lim YY{a,,^)= lim (d,,.) (10) 

j=l k=l 

which (for given coupling constants) is a function of the magnetic field H 
and the temperature T. The second equality in Eq. (10) follows from trans- 
lational invariance, which guarantees that the local magnetization m- j = 
limjv,M->oo('7j ^) is the same for all lattice points, while (cj^) still depends 
on the system size. The magnetization can also be obtained as a derivative 
of the free energy f{H,T) by m{H,T) = -df{H,T)/dH. For sufficiently 
high temperature {T > T^), m{H,T) is a smooth monotonously increasing 
function of the field, but at a certain temperature (the critical ov Curie 
temperature), the derivative dm{H,T^) /dH (the magnetic susceptibility) 
becomes infinite at H = 0, and for T < the function m{H,T) has a 
jump discontinuity at H = 0, see figure 1.1 in (Baxter, 1982). The jump 
corresponds to a spontaneous magnetization 

mJT)= lim m(H,T) = - lim m(H,T) . (11) 

' H^o_^ ^ ' H^o_ ^ ' ' ^ ' 

The critical temperature is determined by 

sinh(2K) sinh(2L) = 1 , H = , (12) 

which is just the condition that the Ising model is self-dual, i.e., invariant 
under duality transformation (Kramers and Wannier, 1941). 

The spontaneous magnetization of the system vanishes in the disordered 
phase (r > T^) and is finite in the ordered phase {T < T^), behaving as 
mQ^T) ~ [T^ — T)^ with the critical exponent^ [i = 1/8 as the critical point 
is approached [T < T^), see Figure 4. It shows also the specific heat 

c,iT) = -t'-^I^ (13) 

at constant magnetic field H, which displays a logarithmic divergence at the 
critical point (Onsager, 1944). Thus the corresponding critical exponent is 
a = 0. Other critical exponents can be defined for the singular behaviour 



^We stick to conventional notation since confusion between critical exponent /3 and 
inverse temperature /3 is unlikely. 
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Figure 4- Spontaneous magnetization and specific heat of tlie isotropic Ising model 

of physical quantities near critical points, see (Baxter, 1982) for details. 
2.3. ISING QUANTUM CHAIN 

The Ising quantum chain is defined by the Hamiltonian 

. / N N \ 

^iV = -2 lE^.-J-J+i + E-ll (14) 

acting on the tensor product space 0jLi — C^^. Here, cij is defined as 

(7^ = I® I® . . .0 I®(7^ ® I® . . .0 I® I (15) 

J ^ ^ ^ ^ ^ ^ ^ ' 

j — 1 factors N — j factors 

and (7j accordingly, where and denote Pauli matrices and I the iden- 
tity, which in the usual basis have the form 

'-{ID- 

The Hamiltonian (14) can be viewed (see (Fradkin and Susskind, 1978; 
Kogut, 1979) and (Baxter, 1982, pp. 266-267)) as anisotropic limit of the 
transfer matrix of a layered zero-field square lattice Ising model at finite 
temperature, where the "horizontal" coupling term in Eq. (9) is replaced by 
^j'^jk'^j+ik (S'l^d /i = 0), by performing the simultaneous limits Kj — ?■ 
and L — 7- CO while keeping £j = /i'jexp(2L) finite. Therefore, the zero- 
temperature properties (i.e., eigenvalues and eigenvectors) of the quantum 
mechanical system defined by (14) correspond to the finite-temperature 
behaviour of the classical 2D Ising model. In particular, the critical line 
(12) of the square lattice Ising model corresponds to £j = 1 in this limit. 

The Ising quantum chain can be treated by similar methods as the 2D 
Ising model. The most popular technique following Lieb, Schultz and Mattis 
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(1961) employs a Jordan- Wigner transformation to map the model onto a 
free-fermion system, and thereby reduces the problem of diagonalizing the 
Hamiltonian (14) (which is a 2^ X 2^ matrix) to a diagonalization of a 
quadratic form (which amounts to diagonalizing an A^xA^ matrix). From 
the set of its N eigenvalues {fermion energies ov frequencies) , any eigenvalue 
of (14) can be obtained summing the elements in one of the 2^ different 
subsets. 

To be more precise, one has to consider a modified Hamiltonian to avoid 
the appearance of non-local boundary terms in the fermionic language. It 
is given by the direct sum of the even part (under spin reversal) of the 
Hamiltonian (14) and the odd part of the same Hamiltonian with antiperi- 
odic boundary conditions (i.e., cr^_^_i = —<^i)- This is the mixed sector 
Hamiltonian , for details see (Baake et al., 1989). 

2.4. CONFORMAL SYMMETRY 

Two-dimensional lattice systems at criticality are closely related to con- 
formal field theory (Belavin et al., 1984; Cardy, 1987; Ginsparg, 1990; 
Cardy, 1990; Christe and Henkel, 1993), for a more mathematical introduc- 
tion to conformal field theory the reader is referred to (Gawedzki, 1989). 
This implies that the largest eigenvalues Ai(A^) > A2(A^) > ... of the trans- 
fer matrix are asymptotically degenerate at the critical point. The explicit 
behaviour for an isotropic system is (with j > 1) 



Conformal symmetry predicts the scaling exponents Xj through irreducible 
representations of the Virasoro algebra 



Roughly speaking, Lq is related to the logarithmic derivative of the transfer 
matrix (through a proper scaling limit) and the with negative integers m 
act as ladder operators for the scaled gaps. Explicitly, one gets Xj = A-|-A-|- 
r -\-f where A and A are highest weights of two irreducible representations 
and r, f G Nq. The generator c in Eq. (18) commutes with all L„ and is a 
complex number in any irreducible representation. It is called the central 
charge of the theory and plays a crucial role in a better understanding of 
the notion of universality. The value of the central charge is related (Blote 
et al., 1986) to the finite-size corrections of the largest eigenvalue Ai by 




(17) 



[Lm, Ln] = (m - n)Lm+n + —m{m^ - 1)5, 



m,nel.. (18) 



log (Ai) = -/3/iV + — + o(iV-i), 



(19) 
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where / denotes the bulk free energy per lattice site. 

The critical Ising model is described by the simplest non-trivial unitary 
minimal conformal field theory with central charge c = 1/2. In this case, 
there are three irreducible representations which have highest weights A G 
{0, 1/2, 1/16}. The generating function (or character) of the corresponding 
modules (obtained as Xa(-^) ~ tr(2:^'')) can be written (Baake, 1988) as 



OO 

16 

m = l 



= 1 + z^ + z^ + 2z* + 2z^ + . . . 

= Z2 {l + z + z'^ + z^ + 2z* + ...) ^20) 

= zT^ {1 + z + z"^ + 2z^ + 2z'^ + . . .) 



where 11^(2:) is the generating function of the number of partitions 

oo 

nv(-)=nrT^- (21) 

m = l 

The scaling (or universal) part of the Ising model partition function for pe- 
riodic boundary conditions is a quadratic form in characters, and comprises 
the representations with conformal weights (A, A) = (0, 0) (corresponding 
to the identity operator), (A, A) = (1/16, 1/16) (spin-density operator) and 
(A, A) = (1/2,1/2) (energy-density operator). 

For the Ising quantum chain (14), which alternatively can be viewed 
as the logarithmic derivative of the transfer matrix of the square lattice 
Ising model with respect to the spectral parameter, analogous relations to 
Eqs. (17) and (19) hold for the energy gaps Ej — Ei and the finite-size cor- 
rections of the ground-state energy Ei^ respectively. Here, the eigenvalues 
Ej of the Hamiltonian play the part of — log(Aj). In addition, the scaling 
factor of the universal term has to be modified to correct for the anisotropy 
introduced by the Hamiltonian limit. For the mixed-sector Hamiltonian of 
the Ising quantum chain mentioned at the end of the previous subsection, 
the periodic sector contains weights (A, A) = (0,0), (0,1/2), (1/2,0), and 
(1/2, 1/2), which gives the field theory of a free massless Majorana fermion. 

3. Aperiodic Ising Chains 

It is now time to turn to aperiodic Ising models, commencing with the 
simplest case of the classical ID Ising chain (1). Here, we restrict ourselves 
to chains with a constant magnetic field Hj = H and two ferromagnetic 
(i.e., positive) coupling constants Jj G {</a, j = 1,...,A^. The actual 
sequence of coupling constants is determined by a two-letter substitution 
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rule on the alphabet {a,b}. Substitution sequences are treated extensively 
in other parts of this volume (see also (Allouche and Mendes France, 1995)), 
thus we do not need to go into detail here. 

A number of investigations has been performed on this or related mod- 
els, see the references listed in the Bibliography. In particular, exact renor- 
malization group transformations were used to study the ground-state and 
thermodynamic properties of Fibonacci Ising chains without and with a 
uniform magnetic field (Achiam et al., 1986), also for a mixture of ferro- 
and antiferromagnetic couplings (Tsunetsugu and Ueda, 1987). The zero- 
temperature ground-states for the case of uniform couplings and quasiperi- 
odic field have also been investigated (Luck, 1987; Sire, 1993). 

Of course, it is still impossible to observe a phase transition at non-zero 
temperature in these ID systems, but nevertheless interesting properties 
arise. As an example, we consider the Lee-Yang zeros of the partition func- 
tion in the field variable uj = exp{2(3H) and their distribution on the unit 
circle, compare (Baake et al., 1995; Simon et al., 1995) for complementary 
material. 

As one of our main examples throughout this article we choose the silver 
mean? substitution rule g 

^^nL' ''^^HD' AJ = 1±V2. (22, 

Here, is the corresponding substitution matrix with eigenvalues A^. By 
iterated application of g on the initial word Wq = a one constructs words 
'^n+i = e{w„) = w„w„_iw„ of increasing length 

\w„\ = f„ + f„-i , (23) 

where the silver mean numbers are defined through 

/„+i = 2/„ + /„_!, /o = 0, /i = l. (24) 

Two consecutive numbers are coprime, and /„_|_i//„ — ?■ A+ for n — ?■ oo. 

To any finite word w in the alphabet {a, b} we associate an Ising chain 
(with periodic boundary conditions) with N = \w\ sites by choosing the 
coupling constant Jj in Eq. (1) to be (J^) if the jth letter of w is an a 
(5), respectively. Thus, we have two different local transfer matrices 

T, = {iozX'''{ (25) 



^The irrational number = 1 -|- t/2 = [2; 2, 2, 2,2,.. .] is called the silver mean for 
obvious reasons in comparison with the golden number r = {l + \/E)/2 = [1; 1, 1, 1, 1, . . .]. 
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Figure 5. Lee-Yang zeros and their integrated density on the unit circle with gap labels 



where uj = exp{2(3H), Zg = exp (2/3/^) and s G {a,b}. This way, we ob- 
tain periodic approximants for the silver mean Ising chain. Denoting the 
corresponding product of local transfer matrices by T„, they satisfy the 
recurrence (ra > 0) 



T = T T T 

' n-\-l ' n ' n — 1 ' n 



(26) 



with Tq = and Ti = T^,. One can also derive a recursion relation for the 
partition function tr(7^) using the so-called trace map, see (Baake et al., 
1993) for details. Note that the partition function is real for = 1 (and 
real z^) because T* = a^T^a^ with given by Eq. (16). 

From Eq. (25) it is obvious that the partition function is essentially 
a polynomial in the three variables cj, z^ and z^. The quite general Lee- 
Yang theorem (Ruelle, 1969) guarantees (for ferromagnetic couplings, i.e., 
Zg > 1 for s = a,b) that all zeros in the complex variable uj lie on the 
unit circle = 1 corresponding to imaginary values of the magnetic field. 
Figure 5 shows the zeros and their integrated density on the unit circle for 
the periodic approximant Wq (of length |wg| = 99) with couplings z^ = 3/2 
and Zf^ = 100. In contrast to the periodic case (see above), the distribution 
in Figure 5 shows a characteristic gap structure reminiscent of the general 
picture for the electronic problem which is well understood in terms of 
the famous gap labeling theorem (Bellissard, 1986; Bellissard et al., 1992). 
Analogous results have also been observed for the Fibonacci (Baake et al., 
1995) and Thue-Morse (Simon et al., 1995) Ising chains. 

For the finite system attached to w^, where the integrated density on 
the unit circle is approximated by the spectral step function, we obtain 
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plateaus at values 



in 

< m </„ + /„_! , (27) 



fn + 



by construction. These are independent of the particular choice of the cou- 
pling constants. Replacing the finite system of length p (with periodic 
boundary conditions) by an infinite system of period p, the integrated den- 
sity becomes an absolutely continuous function with the plateaus (and the 
values on them) persisting. For length p = 2, this is explicitly shown in 
the Appendix; for the general statement, one employs the idea of Bloch's 
theorem. 

Since /„ and /„_i are coprime, one can write (Baake et al., 1993) m in 
Eq. (27) as 

m = fi-f„ + p-f„_, (28) 

with i^jV E Z (they can be made unique by proper restriction). Taking 
(yU, v) as label, we select a sequence of gaps in successive approximants 
which converge both in size and position. In our example of Fig. 5, the 
labels corresponding to the largest gaps are also given. In the limit ra — ?■ oo 
(which approaches our aperiodic chain) we obtain the dense set of gap 
values 



(29) 



This coincides precisely with Bellissard's result for the electronic spectra of 
silver mean Hamiltonians. In the present context of magnetic chains, the 
meaning is that the derivative of the specific free energy, integrated along 
closed curves in the cj-plane which do not hit the set of zeros, is quantized. 
More precisely, 

/ f'iio) duj (30) 



27ri Jc 

is a number of the set (29) if C is a path that starts at the origin, and 
returns to it through a gap on the unit circle. 

Though we have sketched which gaps may appear, one still has to prove 
that they do not close in the limit. This can probably be shown with a 
slightly modified version of the constructive gap labeling of (Raymond, 
1995) combined with the approach of the Appendix. 



4. Aperiodic Ising Quantum Chains 

In contrast to the ID classical chain, the 2D Ising model has a phase tran- 
sition at finite temperature. Consequently, as follows from the discussion 
of Section 2.3, the Ising quantum chain (14) inherits this phase transition 
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with the temperature of the 2D model replaced by an expression in the 
coupling constants. 

The Ising quantum chain and related models (notably the XY quan- 
tum chain (Luck and Nieuwenhuizen, 1986; Satija and Doria, 1988)) with 
an aperiodic sequence of coupling constants have been studied frequently, 
a list of articles is given in the Bibliography. It also contains a number of 
papers devoted to the related problem of a 2D Ising model with layered ape- 
riodicity. For both cases, most of the early results concentrate on the special 
examples of the Fibonacci (Doria and Satija, 1988; Igloi, 1988; Tracy, 1988; 
Benza, 1989; Ceccatto, 1989; Henkel and Patkos, 1992) and the Thue-Morse 
(Doria et al., 1989; Lin and Tao, 1990) Ising chains which very closely resem- 
ble the periodic situation. However, on the basis of two different examples 
of three-letter substitution rules, Tracy (1988) already conjectured that 
the critical properties coincide with those of the periodic model only if the 
corresponding substitution matrix has just one single eigenvalue of mod- 
ulus greater than one. This has later been substantiated by investigation 
of models built on generalized Fibonacci sequences (Benza et al., 1990; 
You and Yang, 1990), where also the influence of particular choices of 
initial conditions of the substitution were noticed (Lin and Tao, 1992; 
You et al., 1992). 

Finally, it was shown by scaling and perturbative arguments (Luck, 
1993, 1994; Grimm and Baake, 1994) that the the critical behaviour de- 
pends on the fluctuations in the sequence of coupling constants, resulting in 
the same behaviour as the periodic case only if the fluctuations are bounded. 
In particular, this is the case for substitution rules for which all but the 
largest eigenvalue of the substitution matrix lie inside the unit circle. If the 
characteristic polynomial is irreducible over the integers, this is equivalent 
to the statement that the largest eigenvalue of the substitution matrix (the 
Perron-Frobenius eigenvalue, which is assumed to be larger than one) is a 
Pisot-Vijayaraghavan number or PV number, for short. These are real al- 
gebraic integers > 1 with the property that all their algebraic conjugates 
(except 1? itself) lie inside the unit circle. 

As an example, we consider sequences defined by the substitution rules 

Sk--l^l^k, Mk={l\), A± = l(l±v^fcTl) . (31) 

For A; = 1, this is nothing but the familiar Fibonacci sequence and = r = 
(l-|--v/5)/2 = [1; 1, 1, 1, 1, • • •] is the golden mean which is a PV number. The 
second case (k = 2) is marginal, with eigenvalues = 2 and A^ = — 1, 
whereas for A; > 3 both eigenvalues |A^| > 1. As above, we define our 
sequences by iterated applications of the substitution rule w^_^^ = Qj^{w^) 
on the initial word = a. 
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Figure 6. Fluctuations of letter frequencies in substitution sequences (31) for = 1, 2, 3 

The second eigenvalue governs the behaviour of fluctuations in the 
sequence. To see this, consider the deviations from the mean frequency 
= 1 — 1/A^ of the letter a 

g{N) = \w^{N)l-p,N , g^=g{\w„\), /i„= max b(iV)|, (32) 

where w^{N) denotes the word obtained by truncating the limit word 
(fixed point of the substitution) after N letters, and \w\^ counts the 
number of occurrences of the letter a in the word w. Generically, \g^\ (and 

hence /i„ > |^„|) grows as the rath power of \ In the marginal case, \g^\ 
is bounded, but /i„ may still diverge logarithmically with the length of the 
word (i.e., linearly in ra). In the PV case, we have hounded fluctuations^ 
i.e., \g^\ —7- for ra — ?■ oo and /i„ is bounded (Bellissard et al., 1989). In 
Figure 6, this behaviour can clearly be seen for the substitution sequences 
(31) with k = 1,2,3. 

If we consider the Ising quantum chain (14) with two different coupling 
constants and arranged according to the sequence of letters a and b 
in a substitution sequence, the condition for criticality reads (p^ = 1 — p^) 

Pa log(ea) + Pb log(efc) = (33) 

which can be solved by parametrizing the couplings by 

e„ = r-P" , = (34) 

in terms of an arbitrary positive real number r (Benza, 1989). As men- 
tioned before, the Hamiltonian (14) can be mapped onto a free fermion 
problem, which means that all its 2^ eigenvalues are obtained as the sums 
of N fermion energies {fermion frequencies) with coefficients or 1 (there 
are obviously 2^ different choices of the coefficients). Although we can- 
not present an analytic solution, this connection allows for a very efficient 
numerical treatment. 
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In Figure 7, the integrated densities of the critical fermion frequencies 
(normalized by their maximal values) are shown for the three substitution 
rules under study, for coupling constants given by Eq. (34) with r = 3. Here, 
we present the result for quantum chains obtained from ra = 12 (A; = 1), 
ra = 8 (A; = 2), and ra = 7 (A; = 3) iterations of the substitution rule, 
corresponding to chains of length N = 233, N = 171, and N = 217, respec- 
tively. Again, a characteristic gap structure similar to that observed in the 
Lee-Yang zeros of the ID classical chains is apparent, and we suppose that 
the gap labeling theorem describes this situation as well. One important 
difference, however, occurs at the low-energy edge of the fermion spectrum. 
Compatibility with conformal invariance requires here a linear integrated 
density in the thermodynamic limit. Though this looks plausible in the case 
k = 1, the other integrated densities suffer from the fluctuations. To see 
this more clearly, one has to consider the appropriately scaled low-energy 
spectrum of the quantum chain. 

In the second half of Figure 7, we show the low-lying part of the com- 
plete energy spectrum (of the mixed sector Hamiltonian mentioned above 
in Section 2.3), which has been scaled such that the third gap (i.e., the 
gap between the third excitation and the ground-state energy) has width 
one. For the Fibonacci case {k = 1), one clearly recognizes the equidistant 
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scaled spectrum predicted by conformal invariance for the periodic system; 
a closer examination reveals that also the degrees of degeneracy (in the 
thermodynamic limit) match the prediction as described in Section 2.4. As 
an additional ingredient, one needs a modified scaling factor due to the 
anisotropy of the model. Here, this factor is known analytically and has the 
form N/ [2kv) where v is the fermion velocity v = 2 log(r) /(r — r~^) (Luck, 
1993) with r given in Eq. (34). 

In the other two cases, the spectrum behaves quite differently from 
that of the Fibonacci chain. The normalization needed to keep the scaled 
gap finite grows faster than linearly in the system size N; for A; = 2 it 
contains an additional logarithmic contribution in the system size, for k = 
3 it grows as a power of N. This is in agreement with Luck's criterion 
(1993) for the relevance of aperiodicity — if it has bounded fluctuations, 
the aperiodicity is irrelevant and the critical behaviour is precisely the same 
as in the periodic case. However, if the fluctuations are unbounded, then 
the critical behaviour is indeed affected by the aperiodicity. The scaling 
behaviour in the marginal case with logarithmically diverging fluctuations 
{k = 2) has been investigated in (Berche et al., 1995) for the 2D model 
with layered aperiodicity. 

5. Ising Model on Labyrinths 

The most interesting case where one might hope to obtain analytic re- 
sults is that of a 2D Ising model with an aperiodic distribution of coupling 
constants. As discussed above, the case of layered aperiodicity is covered 
by the Ising quantum chain, and numerous investigations (based on se- 
ries expansion or numerical techniques such as Monte Carlo simulations) of 
Ising models on aperiodic graphs (mostly the Penrose tiling) have been per- 
formed, see the corresponding list in the Bibliography section. In contrast 
to the layered case, most results indicate that the critical behaviour remains 
unaffected by aperiodicity (Godreche et al., 1986; Wilson and Vause, 1988; 
S0rensen et al., 1991). There are, however, specific influences of the local 
order on the Curie temperature (Bhattacharjee, 1994; Ledue et al., 1995) 
and on frustration in antiferromagnetic systems (Okabe and Niizeki, 1988; 
Oitmaa et al., 1990; Duneau et al., 1993; Ledue et al., 1993) while several 
hints on non-universal behaviour of Ising models on quasicrystals (Bhat- 
tacharjee et al., 1987; Minami and Suzuki, 1994) seem non-conclusive and 
need further investigation. First results for 3D Ising models (Bose, 1987; 
Okabe and Niizeki, 1990) show no signs of deviant critical properties due 
to aperiodicity. 

The situation is similar for another class of critical phenomena compris- 
ing percolation of walks on quasiperiodic graphs. A number of publications 
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devoted to this subject are listed in the Bibliography, see in particular 
the article by Briggs (1993) which contains numerical estimates for critical 
points and exponents of self-avoiding walks on several quasiperiodic graphs 
obtained by a series expansion analysis. Since walks are only sensitive to 
the topology of the graph, it is even harder to detect non-universal be- 
haviour. This requires further investigations which go beyond the scope of 
this article. 

Let us return to classical Ising models on 2D graphs. Apart from nu- 
merical approaches, the concept of "Z-invariance" (Baxter, 1978) has been 
applied to obtain analytic results for the eight-vertex and certain IRF mod- 
els (Korepin, 1986, 1987; Choy, 1987, 1988; Antonov and Korepin, 1988) 
on (rhombic) duals of regular de Bruijn grids (de Bruijn, 1981). Here, we 
want to discuss the Ising model on 2D quasiperiodic graphs which are con- 
structed from a rectangular grid. For definiteness, we concentrate on the 
so-called Labyrinth tiling (Sire et al., 1989) for which the underlying grid 
is built on silver mean sequences, compare (Baake et al., 1994). However, 
most results can easily be extended to more general classes of graphs. 

5.1. THE LABYRINTH TILING 

The silver mean chain is obtained by repeated application of the two-letter 
substitution rule g of Eq. (22) to the letter Wq = a, i.e., = q{w^). From 

this, the Labyrinth tiling can be constructed by considering an orthogonal 
Cartesian product of two identical silver mean chains in proper geometric 
representation, i.e., the letters b and a are represented by long respectively 
short intervals with length ratio given by the silver mean A := A+ = 1 + ^2. 
In this way, one arrives at a rectangular grid, and the Labyrinth tiling is 
obtained by choosing one of the two subgrids (defined as the set of vertices 
which can be reached from a chosen vertex by an even number of steps on 
the grid, or, in other words, as the two equivalence classes with respect to 
the corresponding equivalence relation) and connecting neighbouring points 
of the subgrid, see Figure 8. The graph defined that way is topologically 
equivalent to a square lattice. 

Furthermore, we can naturally define a dual graph by connecting the 
points which lie on the other subgrid. Since the words obtained from 
the silver mean substitution rule (22) consist of an odd number of letters 
and are palindromic, the two dual graphs are equivalent to each other (even 
for a finite patch constructed from a word w^) by reflection respectively by 
rotation through 90 degrees (which coincide due to the reflection symmetry 
of the graph with respect to the diagonal). This also implies a one-to-one 
map between the edges (bonds) of the graph and the dual graph; in fact, one 
has two choices due to the reflection symmetry of the Labyrinth. Either, 
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Figure 8. Labyrinth tiling with underlying grid and rapidity lines 

one associates pairs of edges related by rotation through 90 degrees, or 
pairs that are mapped onto each other by reflection at a grid line. Here, we 
opt for the latter version because in this case the intersections between the 
Labyrinth and its dual (see Figure 9) occur between pairs of dual bonds, 
which is exactly what we need for the Ising model duality transformation 
discussed in Section 5.2 below. 

Figure 9 shows the three different types of mutually dual tiles (dashed) 
and vertex configurations that appear in the Labyrinth. Their frequencies 
(in the infinite tiling) can easily be calculated from a induced substitution 
rule, which gives 5 — 2A for the square, 6A — 14 for the kite and 10 — 4A 
for the trapezoid, respectively. Since the tiling has fourfold rotational sym- 
metry (Sire et al., 1989), the different orientations of the tiles (one for the 
square and four for the kite and the trapezoid) and corresponding vertex 
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Figure 9. The mutually dual tiles and vertices of the Labyrinth 

configurations occur equally often. The same is true of the three different 
edge types — long (denoted by i), medium (m) and short (s) — of the 
Labyrinth, their frequencies (normalized relative to the number of vertices, 
i.e., i^i + + /Us = 2 as there are twice as many bonds as vertices) are 
calculated to give /^^ = 1, = 2A — 4 and ji^ = b — 2A, respectively. Here, 
the long and short bonds occur in two, the bonds with medium length in 
four different orientations. Note that an analogous construction proceeding 
from a general two-letter substitution rule has four tile and vertex config- 
urations (but still only three different edge lengths), one is missing here 
because aa is not contained in the silver mean sequence as a subword. 

Now, we want to consider periodic approximants constructed from the 
finite patches described above. This is done using a periodic rectangu- 
lar grid and performing the same construction as before. Since \w^\ = 
2|w„_i| -|- |w7„_2| is odd, we cannot use the grid obtained from repeating 
periodically, because we need the existence of two distinct subgrids. 
Therefore, we instead identify the first and last letter in \w^\ (which both 
are 5, ra > 0), or, in other words, construct our periodic grid by repeating 
the word w'^ {n > 0) periodically (for both horizontal and vertical direc- 
tion), where w'^ denotes the word that is obtained by deleting the last letter 
in w^. In this way, one arrives at a series of periodic approximants with 
periods \w^\ — 1 which first of all do not contain any mismatches and fur- 
thermore are still equivalent to their duals constructed from the other 
subgrid. This follows from the obvious fact that the two periodic sequences 
built on the word w'^ and on its reverse are equivalent by translation. 

5.2. ISING MODEL AND DUALITY TRANSFORMATION 

Given a periodic approximant C^, we define our Ising model by placing 
the Ising variables on the vertices and associating coupling constants to 
the edges connecting them. Taking into account different orientations, we 
have eight different types of bonds, and therefore (in this setup) up to 
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Figure 10. The eight different bonds and their associated coupUng constants 

eight different coupling constants at our disposal. Let us, with [i = l/ksT, 
use K^y/(3 and L^y/(3 {xy G {aa, ab,ba,bb}) to denote the (ferromagnetic) 
coupling constants for bonds which are "raising" and "lowering" diagonals 
of a rectangle with abscissa x and ordinate y in the underlying grid, re- 
spectively, see Figure 10. In what follows, we consider the ferromagnetic 
zero-field case only, thus our parameter space is since only the products 
of the coupling constants with the inverse temperature enter in the compu- 
tation of the partition function. The reflection symmetry of the Labyrinth 
mentioned above implies that we can define a transformation on the space 
of coupling constants through 

\ / \ 

(35) 





which leaves the partition function of the Ising model invariant. 
The duality transformation * is an analytic involution 

(36) 

on M^, where K*y G M_,_ and L*^ G M_,_ are defined by 

sinh(2K:^) sinh(2L,^) = sinh(2L:^) sinh(2K,^) = 1 . (37) 

Under this transformation, the partition function is changed only trivially. 
This can be seen by investigating two exact graphical expansions of the 
partition function in terms of polygons on the Labyrinth {high-temperature 
expansion) respectively its dual {low-temperature expansion) which is again 
the Labyrinth, see (Baake et al., 1994) for details. 

The Ising model is called self-dual if /i' „ = K*, and L^„ = L*„, i.e., 

o t> xy xy xy •^y 



where S^y is defined as 



Saa = Sab = Sba = 5*66 = 1 , (38) 

S,y ■.= smh{2K,y)smh{2L,y) . (39) 
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Consider a smooth curve in parameter space along which all coupling con- 
stants are strictly increasing (ranging from to infinity) and which is 
mapped onto itself by duality. It follows from Peierls' argument (Peierls, 
1936; Ellis, 1985) that along such a curve one encounters at least one phase 
transition point. However, phase transition points correspond to parameter 
values where the free energy fails to be analytic, and thus must be mapped 
onto other phase transition points under duality. If one assumes that there 
is only one phase transition point on the curve (which is a natural assump- 
tion for the Ising model), it necessarily has to coincide with the intersection 
of the curve with the self-dual surface (which exists since the involution * 
is analytic). Therefore, we suppose that Eq. (38) defines a four-dimensional 
critical surface in our eight-dimensional parameter space. 

Although this result is much stronger than those obtained for Penrose 
tilings, where duality arguments together with correlation inequalities yield 
lower and upper bounds for the critical temperature (Bhattacharjee et al., 
1987), it is not completely satisfactory. The reason is that, due to Peierls' 
argument, we expect a critical surface of codimension one dissecting the pa- 
rameter space, or possibly an even more complicated scenario with several 
transitions, though this seems rather unlikely (because there appears to be 
no other natural order parameter besides the magnetization in our model, 
as long as we stay in the ferromagnetic regime). Still, the assumption of 
uniqueness of the phase transition is crucial in our argument — and we 
presently do not know of any rigorous argument to exclude the possibility 
of other transitions. 

5.3. COMMUTING TRANSFER MATRICES 

However, we can employ another approach based on commuting transfer 
matrices. As mentioned in Section 2.2, the Ising model on the square lattice 
can be reformulated as an IRF model, and the Boltzmann weights of the 
IRF model can be parametrized (in terms of elliptic functions of the spectral 
parameter) such that they satisfy the YBE, which implies commutativity of 
the corresponding row transfer matrices (for different values of the spectral 
parameter). Apart from the spectral parameter, the Boltzmann weights 
depend on one further parameter (the elliptic modulus) which basically 
plays the role of a temperature. 

The values of the spectral parameter entering in the Boltzmann weights 
can be viewed as the differences of two so-called rapidities'^ associated to two 
sets of orthogonal lines {rapidity lines) which intersect in the corresponding 

^This terminology originates in scattering theory in 1 -|- 1 dimensions, where the YBE 
impUes factorization of S'-matrices and rapidities are used to parametrize the relativistic 
energy-momentum relation. 
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face of the IRF model (respectively on the corresponding bonds of the Ising 
model). In fact, the values of these rapidities may be varied from one line 
to the other without destroying commutativity of row transfer matrices, 
because commutativity is a consequence of the local Yang-Baxter relation. 
This is the idea behind "Z-invariance" as introduced by Baxter (1978) for 
the eight-vertex model and later applied to a checkerboard Ising model 
(1986). The term "Z-invariance" stems from the fact that the YBE allows 
to "shift" the rapidity lines around without changing the partition function 
of the model (which is customarily denoted by the letter Z). 

To make this more concrete, let us have a look at the example at hand. 
Besides the bonds of the Labyrinth and the underlying grid. Figure 8 also 
shows the rapidity lines (dashed) which form a dual rectangular grid. The 
horizontal (vertical) lines are labeled by variables and (m^ and m^) 
according to the type of edges of the underlying grid which they intersect or- 
thogonally. Note that each pair of one horizontal and one vertical line inter- 
sects in one face of the grid respectively on one bond of the Labyrinth. Now, 
for arbitrary (in general complex) values of the four rapidities u^, m^, u'^, m^, 
we can use the same parametrization for the Ising model coupling constants 
on the Labyrinth (respectively the corresponding Boltzmann weights for the 
IRF model on the underlying grid) as in the periodic or the checkerboard 
case (Baxter, 1986). Since only the differences of rapidities enter in the 
weights, and since we have the elliptic modulus as an additional parameter 
(which is not allowed to vary from one face to another), this yields a four- 
dimensional subspace of our eight-dimensional space of coupling constants 
where the model is solvable (in the sense of commuting transfer matrices). 

The free energy in the thermodynamic limit can be calculated explicitly, 
it is given by the sum of free energies of the four periodic arrangements of 
the rapidities, properly weighted by their frequencies. Also, the local spon- 
taneous magnetization can be computed, it is independently of the position 
in the lattice given by 



where Q is the elliptic modulus entering in the definition of the coupling 
constants. The local YBE finally reduces the calculation of the magne- 
tization to that of the periodic case, where one can use Onsager's result. 
Unfortunately, the more general situation does not provide a new and inde- 
pendent method to calculate the magnetization. In conclusion, the critical 
points on the solvable surface are determined hy = 1, and the criti- 
cal behaviour is exactly the same as in the periodic case, with a magnetic 
exponent (3 = 1/8. 




(40) 
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The actual expression for the parametrization of the coupling constants 
is neither essential nor very illuminating, and the interested reader is re- 
ferred to (Baxter, 1986; Baake et al., 1994) for the technical details of the 
calculation. However, to give an impression of the kind of constraints im- 
posed by integrability, let us have a look at the corresponding conditions 
on the coupling constants. To this end, we introduce a third involution ?, 

(41) 

on the space of coupling constants, where K^y ^ ^-i- ^^'^ ^xy ^ ^-i- ^''^^ 
defined through 

sinh {2K^y) sinh {2K^y) = sinh (2L,y ) sinh (2L,y ) = 1 . (42) 

Then the conditions on the couplings constants can neatly be summarized 
in the requirement that the equation 

^2 _ fr sinh(2Mj) cosh(Mi + M2 + M3 + M4 - 2Mj) 
jj-^ smh{2Mj) cosh{Mi + M2 + M3 + M4 - 2Mj) 

has to be satisfied for all vertices of the Labyrinth, where Q is an arbi- 
trary constant (coinciding with the elliptic modulus used in the explicit 
parametrization) and where Mj {j = 1,2,3,4) denote the four couplings 
around the particular vertex. Note that the value of Q in Eq. (43) has to 
be the same for all vertices. For the Labyrinth, this results in altogether 
five relations between Q and the eight couplings K^y, L^y. As can easily be 
verified, four of the five conditions can be simplified to 

Saa = Sab = Sba = Sbb = ^ (44) 

the fifth, more complicated relation stems from the vertex with one long, 
one short and two medium length bonds. Therefore, the critical subspace 
Q = 1 on the solvable surface is in turn a subspace (of codimension one) of 
the self-dual surface defined by Eq. (38). 

Eq. (43) clearly reveals the drastic constraints which integrability im- 
poses on the coupling constants, and which is particularly apparent in the 
result for the local magnetization (40). For a general set of couplings, one 
expects that the magnetization at a particular vertex depends strongly on 
its local neighbourhood rather than giving the same result for all vertices 
of the graph. This behaviour is due to the fact that the value of Q is the 
same for all vertices, which for instance makes it impossible to make one 
coupling around a vertex weaker without having to compensate this by en- 
larging one of the other coupling constants. This shows that the solvable 
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cases are certainly not representative in general, and that one has to be 
careful if one wants to draw general conclusions from the analysis of the 
solvable subspace alone. 

It should be clear from the discussion above that the method of com- 
muting transfer matrices extends to arbitrary graphs which are constructed 
from a rectangular grid in an analogous way. Moreover, the result is always 
the same — all that matters are the relative frequencies of the different cou- 
pling constants. The critical behaviour remains the same as for the periodic 
Ising model, irrespective of the properties of the underlying sequence. 

5.4. LOCAL MAGNETIZATION 

From the last remark, it is clear that we need to gain more information 
about the critical behaviour of the generic case in order to understand 
the model completely. Unfortunately, one cannot in general compute the 
partition function analytically, and therefore we use numerical simulations 
as a first approach. The simplest quantity one can measure is the local 
magnetization at a particular vertex, and we expect to see that it will in 
general depend on the local neighbourhood of the vertex. 

Here, we restrict ourselves to the case of three different coupling con- 
stants associated to the three different lengths of bonds, i.e., K^^ = L^^ =: 

f^Js, Kab = ^afc = = ha=- f^^m, and K^^ = hb =■ P'h^ ^^cre the 
subscripts s, m, and I refer to short, medium, and long bonds, respectively. 
From Eq. (44), one now finds that the only solvable case among those is 
the periodic Ising model where all the coupling constants are equal. 

For our numerical simulation, we employ the Swendsen-Wang Monte- 
Carlo algorithm (Binney et al., 1992) for the periodic approximant con- 
structed from the word of length \w^\ = 41 which is obtained after five 
applications of the silver mean substitution rule q of Eq. (22) to the initial 
word Wq = a. For this patch of {\w^\ — 1)^ = 1600 sites, we estimated the 
local magnetization at three different vertices, where we chose represen- 
tatives of the three different types of vertex configurations, see Figure 9. 
The result is presented in Figure 11, where the local magnetization at the 
three vertices is plotted as a function of T/T^. Figure 11(a) corresponds 
to the periodic case Jg = Jm = Ji, in Figure 11(b) the ratios of coupling 
constants are given by Js/Jm = 6/5 and Ji/Jm =4/5, and Figure 11(c) 
finally displays the results for Js/Jm = 7/5 and Ji/Jm = 3/5. Here, the 
critical temperature was kept approximately constant by adjusting the 
coupling constants such that their average (per bond) is the same in all 
three cases, which is what we used to normalize the abscissa in Figure 11. 

Clearly, the results show that the local magnetization becomes position- 
dependent once we leave the periodic case. That is exactly what one expects 
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Figure 11. Local magnetization at three representative sites for difTerent couplings 

— by choosing the coupling along the long bonds to be smaller than the 
other two coupling constants one diminishes the local magnetization at the 
vertices surrounded by four long edges. On the other hand, Figure 11 does 
not show any indication that there might be more than one transition; the 
local magnetization at the three locations falls down at approximately the 
same value of the temperature. We believe that it is possible to prove the 
previous statement rigorously by applying suitable correlation inequalities 
(Griffiths, 1972); this will be investigated further. 

6. Concluding Remarks 

In this article, examples of the effects of aperiodicity on mathematical and 
physical properties of one- and two-dimensional Ising models have been 
discussed. Here, the aperiodicity we introduced is by no means random; 
in fact, the sequences we used to define our models (and the thermody- 
namic limit) are obtained from substitution rules and hence are completely 
deterministic. 

One common feature that is apparent is the characteristic gap structure 
that shows up in various quantities for essentially one-dimensional situa- 
tions. Obviously, the observed locations of the gaps obey the general gap 
labeling theorem for Schrodinger operators. Though we cannot present a 
rigorous proof for this statement at present, we believe that a constructive 
approach along the lines of (Raymond, 1995) is possible and hope to report 
on it soon. 

As we see in the example of the Ising quantum chain, the critical be- 
haviour of statistical systems may be changed by introducing aperiodicity. 
In this case, the fluctuations of the sequence of coupling constants deter- 
mine whether the critical properties are affected or not. Of course, the same 
is true of a 2D Ising model with layered aperiodicity. 

However, the situation appears to be different if one considers the 2D 
Ising model with an "isotropic" aperiodicity, for instance with a non-trivial 
rotational symmetry as in the case of the Labyrinth. Here, the available 



APERIODIC ISING MODELS 



29 



analytic and numerical results suggest that the critical behaviour remains 
unaffected for a large class of systems. 

As stated in the beginning, it was not the purpose of this article to give 
final answers to the questions that arise in this context. However, we are 
confident that some of our observations can be made rigorous and help to 
improve the understanding of the mathematical and physical properties of 
non-periodic statistical systems. 
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A. Magnetic field zeros of periodic Ising chains 

Let us consider the partition function = tr(r") of a periodic model 

with a 2x2 transfer matrix T. Assuming that det(r) ^ 0, we define ip with 
sft(Vp) G [0, tt] by 

. ^ tr(r) , , 

COS = — p== . 45 

2Vdet(r) ^ ' 

By induction, one shows that 

Z„ = 2 [det(r)]"/2 cos(nv3) (46) 

and thus the equation Z„ = is equivalent to cos(ra(/j) = which has n real 
solutions ip = (k — 1/2)tt /n in the interval [0, tt]. For n — ?■ oo, the solutions 
fill the complete interval [0, tt] with uniform density. 

From this, one can easily calculate the density of partition function ze- 
ros in the variable u = exp{2(3H) of the ferromagnetic periodic ID Ising 
chain (i.e., Jj = J > and Hj = H m Eq. (1)) on the unit circle, which is 
done in some standard textbooks on statistical mechanics, see e.g. (Pathria, 
1972, pp. 420-424). However, here we would like to treat the slightly more 
interesting example of an alternating Ising chain with ferromagnetic cou- 
pling constants J2J-1 — '^a ^ ^ ^^"^ '^2j — '^b ^ ^ ^^"-^ uniform magnetic 
field H- = H. Denoting the two corresponding local transfer matrices by 



30 



UWE GRIMM AND MICHAEL BAAKE 



and T^, the partition function of the alternating chain of length N = 2n 
is 



z„ = tr[(r„r,)"] = tr(r;,) 



(47) 



and hence the partition function of a periodic chain of length n with 
transfer matrix T^^ = T^Tj^. Introducing the notation = exp(2/3Jg) > 1, 
z^ = exp(2/3J^) > 1 and 9 = —2i(3H, we find 



2{z,z,y^'co8{e)+2{z,z,r'^' 
-1 -^-1 -^-1, 



tr(r„fe) 

det(r„b) = z^Zf^ + {z^Zf^) ^ - z^z^^ - z^^Zf^ 
and therefore one can solve Eq. (45) for cos(^) to obtain 

COSe= ^(1 - Za^){l - Z^^) COs{lp) - {z^Z^,)'^ 

The N = 2n zeros of the partition function are given by 

,(2A;-ll7r. , 



arccos 



(l-^-^)(l-^^-^) cos(- 



2n 



(48) 
(49) 



(50) 



(51) 



n. 



{0 < Of, < it) and by their conjugate zeros 27r — where k = 1,2, 
For ra — 7- CO, they form (if z^ ^ z^ and [i is finite) two separated intervals 
with gap edges 9^ determined by inserting (/j = and (/j = tt in Eq. (50). 
This yields 



arccos 



^{1-Za^){l-Z^^)-{Z,Z^)-' 



(52) 



(0 < < ^~ < tt) and the corresponding conjugates 2^ — 6^ . In contrast 
to the situation in the periodic Ising chain, the partition function zeros of 
the alternating chains have two gaps, one around ^ = and the other in 
the vicinity of ^ = tt. The latter closes if z^ and z^ become equal, the gap 
at ^ = (which is related to the absence of a phase transition) only for 
[i —7- CO, i.e., at zero temperature, where the distribution of zeros becomes 
uniform. 

The integrated density of zeros V{9/2tt) is defined as the limit of the 
ratio of the number of zeros Oj, < 9 and the total number of zeros, which 
is 2n. Since the distribution of the values of Lp is uniform on [0,7r], this 
function has the following form 



V{e/2'K) 



( 



2n 

1 

2 



arccos 



cos(e)+(^„^j-i 



1-V{1- e/2Tr) 



<9 <9+ 
9+<9<9- 

9-<9 <K 
TT <9 <2tt 



(53) 
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and is shown in Figure 2. The corresponding root density is obtained by 
differentiation and diverges at the four gap edges ib^^ as |^ — (ib^^)|^/^ 
{Lee-Yang edge singularity). 

Finally, let us mention that the same technique applies for periodic Ising 
chains with arbitrary period p, as long as the lengths of the chains N = np 
are integer multiples of p. However, the right hand side of Eq. (45) is a 
polynomial of degree p in cos(^/2) which generally makes it impossible to 
solve Eq. (45) for 9 explicitly. 
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